Numerical delicacies associated with the use of isoneutral mixing operators in ocean models

Ocean models usually rely on a tracer mixing operator which diffuses along isoneutral directions. This requirement is imposed by the highly adiabatic nature of the oceanic interior, and a numerical simulation needs to respect these small levels of dianeutral mixing to maintain physically realistic results. This is a key issue nowadays in oceanic numerical models (e.g., Hansen et al., 2011 ”Earth’s energy imbalance and implications”). For non-isopycnic models, rotated mixing operators must therefore be used. In continuation of the pioneering work of Griffies et al. (1998) and Beckers et al. (2000) in this field, we, first, exhaustively present the performance of various space-time discretizations in terms of stability, accuracy, tracer variance dissipation and min-max violations for the harmonic and biharmonic rotated operators. From this study, we can anticipate the flaws of the different schemes in practical situations. Because global climate models are now targeting increasingly higher horizontal resolution, the question of the viability of an isoneutral biharmonic operator is not only relevant for the regional modeling community but also for the ocean climate community. A new way of handling the temporal discretization of this type of operator is thus introduced. This scheme requires only the resolution of a simple one-dimensional tridiagonal system in